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INTRODUCTION 


Although  the  two-dimensional  cascade  analysis  represents  a  simplified 
version  of  the  actual  three-dimensional  flow  field  which  includes  end  wall 
effects,  the  two-dimensional  problem  gives  significant  insight  into  the 
cascade  flow  field  and  obviously  is  a  necessary  first  step  in  developing  a 
three-dimensional  analysis.  Hence,  cascade  analyses  of  various  types  have 
been  a  subject  of  high  interest  in  recent  years.  Among  the  analyses  being 
pursued  are  inviscid  analyses  (e.g.  Refs.  I  and  2),  inviscid  analyses  with 
boundary  layer  corrections  (e.g.  Ref.  3)  and  Navier-Stokes  analyses  (e.g. 

Ref.  4  and  5).  Each  of  these  approaches  are  viable  under  certain 
circumstances.  Inviscid  analyses  can  give  good  predictions  of  the  blade 
pressure  distribution  for  conditions  where  the  effect  of  viscous  phenomena 
upon  the  blade  pressure  distribution  remains  small.  However,  inviscid 
analyses  require  some  method  of  assuming  airfoil  circulation  to  obtain  a 
unique  flow  solution.  In  general  this  specification  is  straight-forward  for 
sharp  trailing  edged  blades  in  steady  flow  where  the  Kut t a-Joukowski 
condition  serves  to  specify  circulation.  However,  in  cases  where  the 
trailing  edge  is  rounded  or  the  flow  is  unsteady  or  trailing  edge  separation 
occurs,  specification  of  a  proper  condition  is  ambiguous.  In  addition  since 
inviscid  methods  are  devoid  of  viscous  effects  by  definition,  they  obviously 
cannot  give  either  heat  transfer  or  viscous  loss  information. 

If  an  inviscid  analysis  is  combined  with  a  boundary  layer  analysis  in 
either  a  strong  interaction  or  weak  interaction  mode,  some  of  these 
limitations  may  be  relieved.  In  cases  where  the  viscous  displacement  effect 
has  an  insignificant  or  only  small  effect  on  the  actual  blade  pressure 
distribution,  an  inviscid  calculation  can  be  used  to  obtain  the  pressure 
distribution  and  a  boundary  layer  calculation  then  made  to  obtain  heat 
transfer  and  loss  effects.  However,  in  many  cases  such  as  when  boundary 
layers  thicken  significantly  or  separate,  the  viscous  displacement  effect  may 
alter  the  pressure  distribution.  This  may  be  a  particular  problem  in 
transonic  flow  where  the  local  pressure  distribution  and  shock  location 
become  very  sensitive  to  small  change  in  the  effective  passage  area.  In 
these  cases  a  strong  interaction  solution  is  required  to  account  for  the 
mutual  effects  of  the  viscous  boundary  layer  and  the  nominally  inviscid  core 
flow. 
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A  sfroii};  interaction  analysis  may  take  the  form  of  either  a  forward 
marching  procedure  or  a  global  iteration.  For  regions  where  the  outer 
nominally  inviscid  flow  is  supersonic  (and  thus  described  by  hyperbolic 
equations)  a  solution  can  be  spatially  forward  marched  in  the  nominally 
streamwise  direction  with  the  inviscid  and  viscous  regions  coupled  on  a 
station-by-station  basis.  Obviously,  the  analysis  is  limited  by  the 
governing  assumptions  which  include  lack  of  viscous  effects  in  the  outer 
region  and  neglect  of  normal  pressure  gradients  and  streamwise  diffusion  in 
the  inner  region.  The  chief  numerical  difficulty  with  this  approach  is  the 
stiff  nature  of  the  coupled  sets  of  equations  which  is  manifested  in  the 
appearance  of  physically  unrealistic  branching  solutions.  In  regions  where 
the  outer  flow  is  subsonic,  the  outer  flow  equations  are  elliptic  in  nature 
and  in  these  regions  forward  marching  in  the  streamwise  spatial  direction  is 
not  possible  and  consequently  a  series  of  viscous  and  inviscid  calculations 
must  be  performed  in  which  each  corrects  the  other  in  a  global  manner. 

Problems  with  interaction  solutions  become  particularly  severe  in 
transonic  flows  where  both  subsonic  and  supersonic  nominally  inviscid  regions 
are  present  and  where  small  viscous  displacement  effects  may  have  a  major 
effect  on  the  blade  pressure  distribution  and  shock  location.  A  final 
difficulty  with  the  interactive  approach  occurs  when  boundary  layer 
separation  appears.  Here  with  an  imposed  pressure  field,  the  usual  steady 
state  boundary  layer  equations  are  unstable  when  solved  as  an  initial  value 
problem  in  space  in  regions  of  reversed  flow.  However,  the  equations  can  be 
marched  in  space  by  suppressing  the  streamwise  convection  term  in  the 
separated  region  (Ref.  6).  Although  this  approximation  allows  the  solution 
to  be  marched  through  separation,  the  approximation  becomes  progressively 
more  inaccurate  as  the  extent  of  the  separation  zone  or  the  magnitude  of 
either  the  normal  or  the  back  flow  velocities  become  large.  Thus,  calculated 
flow  details  which  may  be  important  (such  as  heat  transfer  at  reattachment) 
may  have  significant  error  when  separation  is  present  and  such  an  interactive 
analysis  is  used.  Other  schemes  have  and  are  being  developed  which  solve  the 
interaction  problem  without  encountering  an  instability  by  either  changing 
the  problem  to  initial  value  in  time  or  iteration  space.  Nonetheless  the 
resulting  solutions  still  retain  the  approximations  of  the  boundary  layer 
equations  and  the  inviscid  flow.  However,  as  with  inviscid  flow  solutions, 
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combined  viscous  and  inviscid  solutions  remain  a  valuable  tool  for  those 
classes  of  cascade  problems  where  the  approximations  adopted  are  valid. 

The  final  procedure  currently  available  is  the  solution  of  full 
ensemble-averaged  Navier-Stokes  equations.  Such  an  analysis  has  been  applied 
to  a  variety  of  cascade  flow  fields  by  Shamroth,  Gibeling  and  McDonald  (Ref. 
7),  Shamroth,  McDonald  and  Briley  (Refs.  5,  8-10)  and  by  Shamroth  and 
McDonald  (Ref.  4).  The  use  of  the  full  Navier-Stokes  equations  for  the 
cascade  problem  allows  use  of  a  single  set  of  equations  for  the  entire  flow 
field  and  thus  removes  the  need  for  an  interaction  analysis  to  couple 
different  equation  descriptions  for  different  flow  regions.  The  analysis 
simultaneously  predicts  both  the  blade  pressure  distribution  and  viscous  and 
heat  transfer  effects. 

To  date  this  Navier-Stokes  analysis  has  been  tested  against  a  variety 
of  sets  of  cascade  experimental  data.  These  include  the  turbine  cascades  of 
Turner  (Ref.  11)  and  Hylton,  Mihelc,  Turner,  Nealy  and  York  (Ref.  12),  the 
compressor  cascades  of  Stephens  and  Hobbs  (Ref.  13)  and  Hobbs,  Wagner, 
Dannenhoffer  and  Dring  (Ref  14)  and  the  compressor  rotor  cascade  of  Dring, 
Joslyn  and  Harden  (Ref.  15).  Since  most  experimental  data  is  limited  to 
surface  pressure  data,  the  comparisons  made  are  primarily  between  calculated 
and  measured  surface  pressure  distributions.  However,  in  Ref.  16  comparisons 
are  made  between  predicted  and  measured  surface  heat  transfer  rates  and  Ref. 
17  gives  velocity  profile  comparisons.  In  general,  the  comparisons  are  quite 
favorable  and  indicate  the  Navier-Stokes  approach  to  be  a  viable  and 
practical  technique  for  predicting  cascade  flow  fields  which  does  not  contain 
the  possible  inadequacies  of  more  approximate  approaches.  Details  of  the 
results  are  given  in  the  cited  references. 

Although  these  results  show  the  Navier-Stokes  approach  to  hold  con¬ 
siderable  promise  as  a  design  tool,  at  the  initiation  of  the  current  effort 
the  code  required  user  expertise  to  obtain  an  accurate,  converged  solution  in 
an  efficient  manner.  The  purpose  of  the  present  effort  was  to  revise  the 
code  so  as  to  allow  a  user  with  little  prior  experience  to  utilize  the  pro¬ 
cedure  in  an  effective  manner.  The  changes  made  to  reach  this  objective  are 
discussed  in  the  present  report.  However,  prior  to  a  discussion  of  these 
changes  a  brief  review  of  the  analysis  is  given. 
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II.  ANALYSIS 


The  present  analysis  is  based  upon  a  solution  of  the  ensemble-averaged 
Navier-Stokes  equations  using  the  linearized  block  implicit  (LBl)  method  of 
Briley  and  McDonald  (Ref.  18).  The  equations  are  solved  in  a  constructive 
coordinate  system  (Ref.  4)  with  density  and  the  Cartesian  velocity  components 
being  taken  as  dependent  variables.  The  application  of  the  LBI  method  to  the 
cascade  flow  field  problem  has  been  discussed  in  some  detail  in  Refs.  4,  and 
7-10.  However,  for  completeness  it  will  be  repeated  here  along  with  a  brief 
discussion  of  the  coordinate  system  and  governing  equations. 

Governing  Equations 

The  present  effort  solves  the  time-dependent  compressible  Navier-Stokes 
equations  to  predict  the  cascade  flow  field.  If  the  computational  spatial 
coordinates  are  £  and  D  where 

£=£(*,y,t)  =  ??(x  »y»t)  r=t  (d 

then  the  continuity  equation,  the  x-component  of  the  momentum  equation  and 
the  y-component  of  the  momentum  equation  are  written  as 
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In  Eqs .  (1-3)  x  and  y  are  Cartesian  coordinates,  t  is  time,  u  and  v  are 

velocity  components,  p  is  density,  p  is  pressure,  and  is  the  stress 

tensor  and  Re  is  the  Reynolds  number. 

The  dependent  variables  chosen  for  the  present  formulation  are  the 
density,  p,  and  the  velocity  components,  u  and  w.  Although  the  code  does 
contain  an  energy  equation  and  calculations  have  been  made  with  an  energy 
equation,  most  calculations  have  been  run  with  the  assumption  T° ,  the 
stagnation  temperature,  equals  a  constant.  With  this  assumption,  the 
pressure  is  related  to  the  velocity  and  density  by 


P  =  p  R 


(u2  +  v2A 

2Cp  J 


(4) 


This  approximation  of  constant  stagnation  temperature  has  been  made  solely  to 
conserve  run  time.  Calculations  made  with  an  energy  equation  (e.g,  Ref.  16) 
have  shown  no  problems  in  converging  and  have  given  surface  pressure  and  heat 
transfer  distributions  in  good  agreement  with  data.  However,  inclusion  of 
the  energy  equation  obviously  requires  additional  computer  run  time. 


Numerical  Procedure 


The  numerical  procedure  used  to  solve  the  governing  equations  is  a 
consistently  split  linearized  block  implicit  (LBI)  scheme  originally 
developed  by  Briley  and  McDonald  (Ref.  18).  A  conceptually  similar  scheme 
has  been  developed  for  two-dimensional  MHD  problems  by  Lindemuth  and  Killeen 
(Ref.  19).  The  procedure  is  discussed  in  detail  in  Refs.  18  and  20.  The 
method  can  be  briefly  outlined  as  follows:  the  governing  equtions  are 
replaced  by  an  implicit  time  difference  approximation,  optionally  a  backward 
difference  or  Crank-Nicolson  scheme.  Terms  involving  nonlinearities  at  the 
implicit  time  level  are  linearized  by  Taylor  expansion  in  time  about  the 
solution  at  the  known  time  level,  and  spatial  difference  approximations  are 
introduced.  The  result  is  a  system  of  multi-dimensional  coupled  (but  linear) 
difference  equations  for  the  dependent  variables  at  the  unknown  or  implicit 
time  level.  To  solve  these  difference  equations,  the  Douglas-Gunn  (Ref.  21) 
procedure  for  generating  alternating-direction  implicit  (ADI)  schemes  as 
perturbations  of  fundamental  implicit  difference  schemes  is  introduced  in  its 
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natural  extrusion  to  sysLoms  of  partial  differential  equations.  This 
technique  lends  to  systems  of  coupled  linear  difference  equtions  having 
narrow  h l ock-hnnded  matrix  structures  which  can  he  solved  ofliciently  by 
standard  b 1 ock-o Li mi nat ion  methods. 

The  method  centers  around  the  use  of  a  formal  Linearization  technique 
adapted  for  the  integration  of  in  it ial -value  problems.  The  linearization 
technique,  which  requires  an  implicit  solution  procedure,  permits  the 
solution  of  coupled  nonlinear  equations  in  one  space  dimension  (to  the 
requisite  degree  of  accuracy)  by  a  one-step  noniterative  scheme.  Since  no 
iteration  is  required  to  compute  the  solution  for  a  single  time  step,  and 
since  only  moderate  effort  is  required  for  solution  of  the  implicit 
difference  equations,  the  method  is  computationally  efficient;  this 
efficiency  is  retained  for  multi-dimensional  problems  by  using  what  might  be 
termed  block  ADI  techniques.  The  method  is  also  economical  in  terms  of 
computer  storage,  in  its  present  form  requiring  only  two  time-levels  of 
storage  for  each  dependent  variable.  Furthermore,  the  block  ADI  technique 
reduces  multi-dimensional  problems  to  sequences  of  calculations  which  are 
one-dimensional  in  the  sense  that  easily-solved  narrow  block-banded  matrices 
associated  with  one-dimensional  rows  of  grid  points  are  produced.  A  more 
detailed  discussion  of  the  solution  procedure  as  discussed  by  Briley,  Buggeln 
and  McDonald  (Ref.  22)  is  given  in  the  Appendix  A. 

Artificial  Dissipation 

Since  the  calculations  of  interest  are  often  at  high  Reynolds  numbers 
typical  of  normal  turbomachinery  applications,  it  is  necessary  to  add 
"artificial  dissipation"  terms  to  suppress  spatial  oscillations  associated 
with  central  spatial  differences  approximations.  This  can  be  done  via  a 
dissipative  spatial  difference  formulation  (e.g.,  one-sided  difference 
approximations  for  first  derivatives)  or  by  explicitly  adding  an  additional 
dissipative  type  term.  For  the  Navier-Stokes  equations,  the  present  authors 
favor  the  latter  approach  since  when  an  additional  term  is  explicitly  added, 
the  physical  approximation  being  made  is  usually  clearer  than  when 
dissipative  mechanisms  are  contained  within  numerical  truncation  errors, 
and  further,  explicit  addition  of  an  artificial  dissipation  term  allows 
greater  control  over  the  amount  of  non-physical  dissipation  being  added. 
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Obviously,  the  most  desirable  technique  would  add  only  enough  dissipative 
mechanism  to  suppress  oscillations  without  deteriorating  solution  accuracy. 
Various  methods  of  adding  artificial  dissipation  were  investigated  in  Ref.  8, 
and  these  were  evaluated  in  the  context  of  a  model  one-dimensional  problem 
containing  a  shock  with  a  known  analytic  solution  (one-dimensional  flow  with 
heat  transfer).  The  methods  which  were  considered  included  second-order 
dissipation,  fourth-order  dissipation  and  pressure  dissipation  techniques. 

As  a  result  of  this  investigation,  it  was  concluded  that  a  second-order 
anisotropic  artificial  dissipation  formulation  suppressed  spatial 
oscillations  without  impacting  adversely  on  accuracy  and  could  be  used  to 
capture  the  nearly  normal  shocks  which  are  expected  in  transonic  cascades 
successfully.  In  this  formulation  the  terms 


d  / 

—  (a  ) 

ax  dx  1 

ay  ^ay  ay  ' 

are  added  to  the  governing  equations  where  <J>  =  u,  v  and  p  for  the  x-momentum, 
y-momentum  and  continuity  equations,  respectively.  The  exponent  n  is  zero 
for  the  continuity  equation  and  unity  for  the  momentum  equations.  The 
dissipation  coefficient  dx  is  determined  as  follows.  The  general  equation 
has  an  x-direction  convective  term  of  the  form  a3<J)/3x  and  an  x-direction 
diffusion  term  of  the  form  3  ( b S <}>/ 3x )/ 3x .  The  diffusive  term  is  expanded 

d(b  d<p/ dx)/dx  -  bd2<f)/dxz  +  db/dx  dcfc/dx  (5) 

and  then  a  local  cell  Reynolds  number  Re^x  is  defined  for  the  x-direction 
by 

ReAx  =  |o  “  db/dx|Ax/b  (6) 

where  b  is  the  total  or  effective  viscosity  including  both  laminar  and 
turbulent  contributions  and  Ax  is  the  grid  spacing.  The  dissipation 
coefficient  dx  is  non-negative  and  is  chosen  as  the  larger  of  zero  and  the 
local  quantity  b  (oxRe/\x- 1 )  .  The  dissipation  parameter  0X  is  a 
specified  constant  and  represents  the  inverse  of  the  cell  Reynolds  number 
below  which  no  artificial  dissipation  is  added.  The  dissipation  coefficient 
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dy  is  evaluated  in  an  analogous  manner  and  is  based  on  the  local  cell 
Reynolds  number  ReAy  anc*  8rid  spacing  Ay  for  the  y-direction  and  the 
specified  parameter  parameter  Oy, 

The  question  arises  as  to  the  values  of  ax  and  Oy  which  should  be 
chosen.  This  was  assessed  both  through  the  model  problem  (Ref.  8),  and 
through  calculations  for  a  Jose  Sanz  compressor  cascade  (Refs.  4,  8  and  9). 
These  results  indicated  that  values  of  a  =  .5  which  corresponds  to  a  cell 
Reynolds  number  2  limitation  would  severely  damp  physical  variations. 

However,  when  a  was  set  in  the  range  .025  °  0.05,  which  correspond  to  a 

cell  Reynolds  number  range  between  40  and  20,  spurious  spatial  oscillations 
were  damped  with  no  significant  change  in  the  calculated  results  as  a  was 
varied  in  this  range.  Further,  as  discussed  in  Refs.  4  and  7-10,  the  results 
obtained  showed  good  agreement  with  data.  This  has  since  been  confirmed  at 
several  other  studies  at  Scientific  Research  Associates  such  as  two-  and 
three-dimensional  transonic  nozzle,  flows  (Ref.  23)  where  a  maximum  acceptable 
value  of  a  =  0.10  has  been  noted  for  most  problems.  In  some  cases  where 
spatial  resolution  may  be  marginal  such  as  at  the  leading  edge  of  a 
relatively  sharp  edged  blade,  it  may  be  necessary  to  increase  a  in  this  local 
area.  However,  a  can  be  decreased  to  0.10  or  below  if  computational  grid 
points  are  added  in  this  region. 

Boundary  Conditions 

The  authors'  experience  in  solving  Navier-Stokes  equations  has  indicated 
the  important  role  boundary  conditions  play  in  obtaining  accurate  solutions 
and  rapid  numerical  convergence.  Improper  specification  and/or 
implementation  of  boundary  conditions  can  lead  to  adverse  stability, 
convergence  and  accuracy  properties  for  the  solution  procedure.  In  regard  to 
specification  of  boundary  conditions,  the  Navier-Stokes  cascade  analysis 
follows  the  suggestion  of  Briley  and  McDonald  (Ref.  24)  which  specifies 
upstream  total  pressure  and  downstream  static  pressure  conditions.  For  the 
cascade  system  shown  in  Fig.  1,  AB  and  CD  are  periodic  boundaries  and 
periodic  conditions  are  set  here. 

Specification  of  upstream  and  downstream  conditions  is  somewhat  more 
difficult.  For  an  isolated  cascade,  boundary  conditions  for  the  differential 
equations  may  be  known  at  both  upstream  infinity  and  downstream  infinity. 


8 


However,  since  computation  economics  argues  for  placing  grid  points  in  the 
vicinity  of  the  cascade  and  minimizing  the  number  of  grid  points  far  from  the 
cascade,  the  upstream  and  downstream  computational  boundaries  should  be  set 
as  close  to  the  cascade  as  is  practical.  However,  when  the  upstream  boundary 
is  placed  close  to  the  cascade,  most  flow  function  conditions  on  the  boundary 
will  not  be  known  since  these  will  have  been  changed  from  values  at  infinity 
by  the  presence  of  the  cascade. 

In  the  present  approach,  the  sugestion  of  Ref.  24  is  followed  which  sets 
total  pressure  on  boundary  BC  (see  Fig.  1).  Unless  boundary  BC  is  very  far 
upstream,  the  flow  velocity  along  this  boundary  will  not  be  equal  to  the 
velocity  at  upstream  infinity  since  some  inviscid  deceleration  will  have 
occurred.  However,  as  long  as  the  boundary  is  upstream  of  the  region  of  any 
significant  viscous  or  shock  phenomena,  the  total  pressure  on  this  boundary 
will  be  equal  to  the  total  pressure  at  upstream  infinity.  Hence,  total 
pressure  is  an  appropriate  boundary  condition  realistically  modeling  the 
desired  flow  condition.  In  addition  to  specifying  upstream  total  pressure, 
it  is  necessary  to  specify  the  inlet  flow  angle.  In  the  present  calcualtion, 
a  value  was  assumed  constant  on  the  upstream  boundary  as  a  specified  angle. 
The  third  condition  set  on  the  upstream  boundary  concerns  the  density  and  a 
zero  density  derivative  at  this  boundary  was  specified  as  a  numerical 
treatment  of  the  boundary.  The  downstream  boundary  was  treated  by  setting  a 
constant  static  pressure  as  a  boundary  condition,  and  by  setting  second 
derivatives  of  both  velocity  components  equal  to  zero  at  this  location.  In 
the  present  application,  a  constant  static  pressure  was  set  at  downstream 
infinity,  and  hence  it  is  assumed  that  the  downstream  boundary  is  located  in 
a  region  where  pressure  is  uniform  although  a  nonuniform  specification  is 
permitted.  The  final  boundary  conditions  to  be  considered  are  the  conditions 
along  the  blade  surface.  Here  no-slip  and  no  through-flow  conditions  were 
applied  leading  to  a  specification  of  zero  velocity  on  the  surface.  An 
inviscid  transverse  momentum  equation  was  applied  on  the  surface  leading  to  a 
boundary  condition  approximation  of  zero  transverse  pressure  gradient  being 
appl ied . 

The  second  item  which  must  be  considered  in  regard  to  boundary 
conditions  is  their  implementation.  Both  the  upstream  and  downstream 
boundaries  have  boundary  conditions  associated  with  them  which  are  nonlinear 
functions  of  the  dependent  variables.  These  are  the  specifications  of  total 
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pressure  on  the  upstream  boundary  and  static  pressure  on  the  downstream 
boundary.  These  nonlinear  boundary  conditions  are  linearized  in  the  same 
manner  as  the  governing  equations,  via  a  Taylor  expansion  of  the  dependent 
variables  in  time,  and  then  solved  implicitly  along  with  the  interior  point 
equations.  Although  points  on  the  periodic  lines  and  the  branch  cut  are 
boundaries  to  the  computational  regime,  they  are  interior  flow  fields  points 
and  must  be  treated  as  such.  The  present  technique  replaces  derivatives  at 
these  points  by  central  differences.  In  addition,  in  regard  to  the  periodic 
lines  the  procedure  inverts  a  matrix  with  strict  periodic  boundary 
conditions;  i.e.,  the  periodic  line  values  are  obtained  from  the  implicit 
solution,  rather  than  from  an  extrapolation  or  averaging  procedure  which  uses 
interior  computational  grid  point  values. 


Turbulence  Models 


The  existing  cascade  analysis  contains  three  possible  turbulence 
models.  These  are  (1)  a  mixing  length  model,  (ii)  a  turbulence 
energy-algebraic  length  scale  model  and  (iii)  a  turbulence  energy-turbulence 
dissipation  model.  In  the  mixing  length  model,  the  turbulent  viscosity  is 
related  to  the  mean  strain  via  a  mixing  length,  £,  such  that 


H-j  =  P  l 


du;  duj  \  <3u j 

dXj  dX|  /  dXj 


1/2 


(7) 


where  pT  is  the  turbulent  viscosity,  p  is  the  density,  i  is  the  mixing 
length,  u i  is  the  i^h  velocity  component  and  x^  is  the  i*1*1  Cartesian 
direction.  Summation  is  implied  for  the  repeated  indices.  The  question  now 
arises  as  to  specification  of  i .  For  the  region  upstream  of  the  trailing 
edge,  the  mixing  length  is  specified  in  the  usual  boundary  layer  manner;  i.e. 

where  k  is  the  von  Karman  constant  and  y+  is  the  dimensionless  normal 
coordinate,  yuT/v.  In  boundary  layer  analysis  &max  is  usually  taken  as 
0.096  where  6  is  the  boundary  layer  thickness  taken  as  the  location  where 
u/ue  =  0.99.  However,  this  definition  of  6  assumes  the  existence  of  an 
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outer  flow  where  the  velocity  ue  is  independent  of  distance  from  the  wall 
at  a  given  streamwise  station;  i.e.,  it  assumes  ue  is  only  a  function  of 
the  streamwise  coordinate-  Although  a  boundary  layer  calculation  will  yield 
solutions  in  which  u  approaches  ue  asymptotically  at  distances  far  from  the 
solid  no-slip  surface,  Navier-S tokes  solutions  for  cascade  flow  fields  do  not 
in  general  predict  a  region  where  u  asymptotes  to  a  constant  value. 
Furthermore,  measurements  of  the  flow  also  show  no  such  region  to  exist  in 
general-  Obviously,  a  proper  choice  of  <5  for  the  Navier-Stokes  cascade 
analysis  is  not  straight  forward.  Calculations  made  in  Refs.  4  and  7-9  have 
set  the  boundary  layer  thickness  by  first  determining  umax,  the  maximum 
streamwise  velocity,  at  a  given  station  and  then  setting  6  via 


(9) 


i.e.,  6  was  taken  as  twice  the  distance  for  which  u/umax  =  k]  .  To  date 
when  k[  has  been  taken  as  0.80,  this  value  has  given  reasonable  results  for 
several  cases-  In  general,  the  boundary  layer  development  including  skin 
friction  and  heat  transfer  is  sensitive  to  the  choice  of  ky  whereas  the 
surface  pressure  distribution  is  relatively  insensitive  to  this  parameter. 

The  model  used  in  the  wake  is  also  a  mixing  length  model  in  which  the 
mixing  length  was  made  proportional  to  the  wake  height,  6,  and  a  linear  of 
growth  of  6  with  distance  was  assumed  based  upon  the  classical  free  jet 
boundary  results  (Ref. 25).  With  the  free  jet  boundary  growth  assumption 


(10) 


where  6pg  and  6gg  are  the  pressure  and  suction  surface  trailing  edge 
boundary  layer  thickness  and  xte  is  the  trailing  edge  location.  The  mixing 
length,  £,  was  taken  as  0.26.  When  using  the  mixing  length  model  the  user 
has  the  option  of  setting  a  transition  location  on  suction  and  pressure 
surfaces.  Most  calculations  run  to  date  have  been  run  with  the  mixing  length 
model. 

In  addition  to  the  mixing  length  model  two  alternate  models  are 
available.  The  first  is  based  upon  the  usual  partial  differential  equation 
for  conservation  of  turbulence  energy  used  in  conjunction  with  an  algebraic 
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length  scale.  The  second  is  a  turbulence  energy-turbulence  dissipation  (k-e) 
model  which  solves  partial  differential  equations  for  both  the  turbulence 
energy  and  turbulence  dissipation.  Converged  calcualtions  have  been  run  with 
both  these  models.  Although,  in  principle,  the  more  sophisticated  turbulence 
models  such  as  the  one-  or  two-equation  models  contained  in  the  present  code 
have  the  potential  for  being  more  accurate,  they  require  specification  of 
model  constants  and  model  functional  forms  (e.g.,  see  Ref.  26).  As  indicated 
in  Ref.  26,  most  of  the  currently  used  models  do  not  appear  to  show  a 
definitive  advantage  over  mixing  length  models  for  general  boundary  layers 
developing  in  arbitrary  pressure  gradients.  The  use  of  such  models  is 
currently  the  focus  of  work  being  conducted  at  SRA. 
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CURRENT  EFFORTS 
Objective 


Although  the  Navier-Stokes  cascade  analysis  has  shown  considerable 
promise,  and  given  favorable  comparisons  with  data,  at  the  initiation  of  the 
present  effort  successful  and  efficient  utilization  of  the  code  required  an 
experienced  user.  Successful  implementation  requires  the  user  to  exercise 
two  codes;  the  first  being  a  coordinate  generation  code  and  the  second  being 
the  Navier-Stokes  code.  The  first  objective  of  the  present  effort  was  to 
automate  these  codes  so  as  to  allow  a  relatively  inexperienced  user  to  (i) 
generate  a  viable  coordinate  system  and  (ii)  to  obtain  accurate  solutions 
with  the  Navier-Stokes  code  in  an  efficient  manner.  The  second  objective  of 
the  effort  was  to  restructure  the  Navier-Stokes  code  so  as  to  obtain 
converged  solutions  within  300  seconds  of  CPU  time  or  less  on  a  modern 
computer  without  resorting  to  code  vector izat ion .  Vector izat ion  should  allow 
a  further  order  of  magnitude  reduction.  As  shall  be  demonstrated,  both 
objectives  have  been  met. 


Coordinate  System 


The  coordinate  system  is  an  important  component  of  the  Navier-Stokes 
analysis.  An  inappropriate  coordinate  system  may  lead  to  difficulty  in 
obtaining  a  converged  solution  and  even  exhibit  physically  unrealistic 
predictions  due  to  geometric  truncation  error.  Therefore,  generation  of  a 
viable  system  is  mandatory.  Any  coordinate  system  used  in  the  analysis 
should  satisfy  conditions  of  (i)  generality,  (ii)  smoothness,  (iii) 
resolvability  and  (iv)  allow  easy  application  of  boundary  conditions. 
Obviously,  a  coordinate  system  must  be  sufficiently  general  to  allow 
application  to  a  wide  class  of  problems  of  interest  if  it  is  to  be 
practical.  The  metric  data  associated  with  the  coordinate  system  must  be 
sufficiently  smooth  so  that  the  variation  from  grid  point  to  grid  point  does 
not  lead  to  a  numerical  solution  dominated  by  metric  coefficient  truncation 
error;  it  should  be  noted  that  this  requirement  differs  from  the  requirement 
of  the  existence  of  a  specified  number  of  transformation  derivatives.  The 
coordinate  system  must  resolve  flow  regions  where  rapid  flow  field  changes 
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occur . 


Finally,  coordinates  should  allow  accurate  implementation  of  boundary 
conditions;  for  the  cascade  this  requires  that  the  metric  coefficients  be 
continuous  across  the  periodic  lines  where  periodic  boundary  conditions  are 
to  be  applied.  The  coordinate  system  presently  used  satisfies  all  these 
requirements . 

In  brief,  the  coordinate  system  consists  of  a  set  of  two  families  of 
curves;  the  £  =  constant  curves  such  as  lines  G'G  or  in  Fig.  1  and  the  r\  = 
constant  curves  such  as  ABCD  or  A'ED'  in  Fig.  i.  Under  the  present  effort, 
the  coordinate  generation  process  has  been  automated  and  the  coordinate 
construction  process  is  as  follows.  The  coordinate  system  is  constructed  by 
first  forming  the  inner  loop  A1 ED1  which  includes  the  blade.  This  is 
followed  by  constructing  an  outer  loop  AGBCD  which  consists  of  periodic  line 
AB  and  CD  and  a  frontal  curve  BC .  Both  the  inner  and  outer  loops  are  then 
represented  by  parametric  curves 


x  =  x(s) ,  y  =  y  ( s) 


(ID 


where  the  parameter  varies  from  zero  to  unity.  The  present  coordinate 
generation  process  utilizes  a  multi-part  transformation  for  the  inner  loop. 
First  x  and  y  are  expressed  as  a  function  of  s1,  the  physical  distance  along 
the  curve.  After  normalizing  s'  such  that  its  range  is  between  zero  and 
unity  a  local  paramet rizAt ion  focuses  upon  the  leading  edge  region  and  places 
s"  =  0.5  at  the  leading  edge  point  where  highest  resolution  is  sought.  This 
location  must  be  specified  by  the  user  and  in  general  will  be  in  the  vicinity 
of  the  front  stagnation  point.  This  initial  transformation  is  done 
automatically  once  the  location  where  maximum  resolution  is  required  is 
specified.  Following  this  first  paramerizat ion ,  s  is  related  to  s"  via  a 
hyperbolic  tangent  parameter izt ion  centered  about  the  leading  edge  point 
where  s"  =  0.5. 

In  this  transformation  for  s"  <  0.5 


s 


I 

21 


In 


I  +  ks" 
I  -  ks" 


(12a) 


and  for  sM  >  0.5 


I  I  +  k  ( I  —  sM  ) 
s  =  l_  2T  ln  I  —  k  ( I  —  s‘* ) 


(12b) 
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where 


k  =  2  tanh 


(12c) 


and  the  parameter  T  controls  the  grid  distribution.  Increasing  T  packs 
points  in  the  vicinity  of  s"  =  0.5 

The  outer  loop  is  then  parameterized  so  as  to  relate  points  on  the  outer 
loop  to  corresponding  points  on  the  inner  loop.  For  two  points  which  are 
periodic,  such  as  G  and  K  in  Fig.  2  to  maintain  periodicity  it  is  necessary 
that  S q  =  1-S^.  This  will  assure  that  coordinate  points  on  the  upper  and 
lower  periodic  lines  will  be  periodically  aligned. 

For  points  downstream  of  the  five  per  cent  chord  location,  the  inner 
loop  and  outer  loop  points  are  made  to  correspond  as  follows.  If  G'  and  K' 
in  Fig.  2  are  points  on  the  inner  loop  aligned  with  G  and  K,  then 


s0=  0.5  SG<  +  o.5(i-sk.) 

V  10 ‘so 


(13) 


Downstream  of  the  five  per  cent  chord  location,  a  locally  polar 
parameterization  is  used  as  shown  in  Fig.  3.  In  this  region  an  origin  for 
the  polar  coordinates  is  chosen  at  xQ  =*  0.05c  and  yQ  halfway  between  the 
upper  and  lower  surface  locations  at  the  five  per  cent  chord  location.  The 
inner  loop  parameterization  s  is  then  tabulated  as  a  function  of  0.  The 
outer  loop  parameterization  then  follows  that  given  Eq .  (13)  to  relate  points 
P  and  Q  on  the  outer  loop  to  P1  and  Q'  on  the  inner  loop.  Finally,  the  outer 
loop  parameterization  is  smoothed  in  the  region  betwen  region  I  (x/c  >  .05) 
and  region  II  (x/c  <  .05).  If  N  pseudo-radial  lines  are  to  be  used,  the  grid 
points  on  each  loop  are  chosen  at  values 


Sj  *  ( i  -  I )  /  ( N-  I )  i  =  0,  N-l 


(14) 
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and  points  having  the  sane  value  of  on  the  inner  and  outer  loops  define 
the  end  points  for  a  given  pseudo-radial  line.  (e.g.  GfG  of  Fig.  2).  Since 
the  loops  are  parameterized  so  that  s  varies  rapidly  in  the  region  of  the 
leading  edge,  this  process  effectively  packs  points  into  the  leading  edge 
region  about  the  location  where  s  =  0.5. 

Following  the  construction  and  parameterization  of  the  inner  and  outer 
loop,  two  intermediate  loops  are  constructed  as  shown  in  Fig.  4.  The  first 
intermediate  loop  is  constructed  downstream  of  x/c  =  0.05  a  given  vertical 
distance,  hi,  from  the  inner  loop.  Upstream  of  x/c  =  0.05,  the  distance  is 
measured  along  a  pseudo-radial  line  as  shown  in  Fig.  4.  A  similar  loop  is 
constructed  inside  the  outer  loop  with  the  vertical  distance  being  maintained 
to  the  end  of  the  periodic  region,  points  R  and  S.  The  inner  secondary  loop 
is  parameterized  to  correspond  to  the  primary  loop  points  it  is  connected 
with  via  the  vertical  or  pseudo-radial  lines;  i.e., 


Similarly, 


's' 

T* 


These  four  parameterized  loops  allow  the  construction  of  the  pseudo-radial 
lines  such  as  G'G  of  Fig.  2  via  the  multi-loop  method  originally  developed  by 
Eiseman  (Ref.  27). 


The  multi-loop  method  requires  introduction  of  a  position  vector  P(r,s) 

with  components  (x,y)  which  will  represent  the  coordinate  location  of  grid 

points.  Based  upon  the  four  loop  construction  process,  vectors  P^(s)  are 

defined  with  i  =  1,2, 3, 4.  Each  P^  has  coordinate  (x£,y-[)  associated 

with  it  at  specific  values  of  s  through  the  previously  described  construction 

process.  A  radial  parameter,  r,  is  then  introduced.  This  parameter  is 

defined  at  the  downstream  boundary,  see  Fig.  4,  as  the  distance  from  the  loop 

in  question  to  the  inner  loop  normalized  by  the  distance  from  the  outer  loop 

to  the  inner  loop.  Thus,  r^  =  0,  r£  =  hi/fi3,  *"3  =  (h  3  -  h2)/h3,  r4  =  1. 

With  the  definition  of  these  quantities  the  general  position  vector  P(r,s)  is 

^ 

related  to  the  loop  position  vectors  Pi(s),  ^(s),  P3(s)  and  P^Cs)  via 


1 


P  (r,s)  =  ( I  —  r)2 ( I  —  a,r )  P, (s)  +  (a,+  2)(l-r)2  rP2(s) 

+  r2[l-a2(l-r)’^(s)  +  (a,  +  2) r2  ( l-r) P3(s) 


05) 
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where 


2 


Q 


a 


2 


3(l-r2)-l 


(16) 


it  should  be  noted  that  at  r  =  0,  P(0,s)  =  PjCs)  and  at  r  =  1,  P(l,s)  = 

-V 

P4  ( s )  .  Further  since  at  r  =  0, 


^  (0,s)  =  [P2(s)-P,(s)](a,+  2) 

and  at  r  =  1 

*  [fys)-P,(s)]«is+2> 


07) 


(18) 


spec i f icat ion  of  the  derivatives  at  the  inner  and  outer  boundaries  is 
determined  by  the  parametric  representation  of  intermediate  loops  2  and  3. 
Thus  the  four  loop  method  allows  specification  of  the  boundary  point 
locations  and  coordinate  angles  at  these  boundaries. 

After  loops  2  and  3  are  parameterized  to  satisfy  the  coordinate  angle  at 
the  boundary  points,  the  grid  is  constructed  as  follows.  If  the  grid  is  to 
contain  M  pseudo-radial  lines  (such  as  line  G'C  of  Fig.  1)  and  N 
pseudo-azimuthal  lines  (such  as  line  JLM) ,  the  values  of  the  pseudo-radial 
coordinate  are 

r(i)  »  i /( N  —  1 )  i«  0,1 ,2,  ...  ,  N-l 

and  the  values  of  the  pseudo-azimuthal  coordinate  are 

S  ( j )  -  j/(M-l)  j  •  0,1,2, ... ,  M-l 

Then  the  position  vector;  i.e.,  the  grid  locations  x,y  for  each  point  in  the 
grid  is  given  by  Eq.  (15). 
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The  preceding  has  assumed  a  uniform  spacing  in  the  radial  direction.  If 
radial  grid  point  concentration  is  desired,  it  is  simply  necessary  to  assume 
a  radial  distribution  function.  The  present  analysis  assumed  a  distribution 
function 


tanhD(l-r) 

tanhD 


(19) 


which  concentrates  points  in  the  wall  region.  Grid  points  are  then  chosen  at 
r(i)  =  (i)/(N-l)  and  the  analysis  proceeded  as  outlined. 

The  input  required  for  the  grid  construction  process  is  as  follows: 


1.  Number  of  pseudo-radial  points 

2.  Number  of  pseudo-azimuthal  points 

3.  Vertical  distance  from  branch  cut  to  periodic  line,  I13 
(See  Fig.  4) 

4.  Angle  of  branch  cut  line,  02  ”  This  is  zero  in  Fig.  4 

5.  Angle  of  leading  edge  line,  O3 

6.  Approximate  leading  edge  radius  of  blade  -  This  sets  the 
pseudo-azimuthal  spacing  in  the  leading  edge 

7.  Reference  flow  Reynolds  number  based  on  free  stream  conditions  and 
blade  chord  -  This  sets  normal  grid  spacing  in  the  vicinity  of  the 
wal  1 

8.  Data  points  or  functional  form  to  define  the  blade  shape 


Grids  generated  with  this  procedure  are  shown  in  Figs.  5  and  6.  Fig.  5 
represents  a  compressor  cascade  corresponding  to  the  cascade  of  Hobbs,  et  al 
(Ref.  14).  This  cascade  was  used  in  the  ensuing  calculations  which  de¬ 
monstrate  the  increase  in  efficiency  of  the  Navier-Stokes  code  obtained  under 
this  effort.  The  second  cascade  which  is  shown  in  Fig.  6  corresponds  to  the 
turbine  cascade  of  Turner  (Ref.  11).  Both  coordinate  systems  were  generated 
from  the  input  given  in  items  1-8  above  and  had  the  required  properties  in 
terms  of  grid  smoothness  and  resolution  to  be  viable  coordinate  systems. 
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The  automated  procedure  developed  to  date  is  confined  to  ,CI  grids  such 
as  those  shown  in  Figs.  6  and  6.  Recently  SRA  has  developed  an  'O'  grid 
generator  (Ref.  16)  and  although  this  procedure  currently  requires 
considerable  user  interaction  to  obtain  a  viable  grid,  automation  of  the  *0* 
grid  gnerator  should  be  straight forwad  based  upon  the  techniques  developed 
under  the  present  effort. 

Navier-Stokes  Analysis 

Having  generated  a  viable  coordinate  system,  the  second  step  in  the 
process  is  to  obtain  a  solution  of  the  Navier-Stokes  equations  for  the 
specified  coordinate  system  and  desired  flow  conditions.  Although  the  deck 
has  proven  capable  of  obtaining  accurate  solutions  for  a  variety  of  cascade 
configurations  and  flow  conditions  (e.g.  Refs.  4,  5,  8,  9,  16),  successful 
operation  had  required  an  experienced  user.  Furthermore,  although  computer 
run  times  were  promising,  they  obviously  could  be  improved,  therefore,  two 
subtasks  were  undertaken  in  the  present  effort.  These  were  (i)  decrease  of 
run  time  and  (ii)  automation  of  input /out put . 

In  regard  to  the  first  item,  decrease  in  run  time,  the  original  casade 
code  was  actually  a  general  code  which  could  be  used  to  solve  a  wide  variety 
of  problems  with  little  if  any  code  revision  required  by  the  user.  Although 
this  generality  represents  a  major  advantage  in  a  general  research  code,  it 
does  exact  a  price  in  terms  of  run  time.  Under  a  previous  effort  some  of  the 
options  originally  available  in  the  general  code  which  were  not  required  for 
the  cascade  problem  were  removed  to  create  a  specific  cascade  version  of  the 
desk.  These  items  included  features  such  as  a  polar  coordinate  option, 
a  two-phase  and  reacting  flow  option,  etc.  Under  the  present  effort  this 
work  was  continued  to  further  decrease  run  time. 

Among  the  items  performed  were  creation  of  specific  tridiagonal  block 
matrix  inversion  routines  for  systems  of  three  equations  and  three  unknowns, 
incorporation  of  small  subroutines  within  the  calling  routines  and  extensive 
recoding  of  frequently  called  subroutines  to  eliminate  small  "DO  LOOPS". 

These  modifications  have  not  changed  the  flexibility  or  generality  of  the 
code;  i.e.,  they  are  transparent  to  the  user  and  have  reduced  run  time  by  a 
factor  slightly  greater  than  two.  In  its  current  state,  the  code  solves  the 
continuity  and  two  momenta  equations  for  turbulent  flow  in  approximately 
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0.00118  secs  per  time  step  per  grid  point  on  a  CRAY-1  computer.  Although  an 
energy  equation  can  be  solved  on  request  by  user,  the  necessary  fast  4x4 
matrix  inverters  have  not  yet  been  written  and,  therefore,  at  present  general 
inversion  routines  must  be  used.  However,  inclusion  of  specific  4x4  matrix 
inverters  is  a  straight-forward  if  somewhat  tedious  task  which  could  easily 
be  part  of  a  future  effort.  It  should  be  noted  that  the  code  is  not 
vectorized  and  careful  vector izat ion  would  be  expected  to  decrease  run  time 
by  a  factor  of  ten.  In  regard  to  run  time  per  case,  successful  calculations 
have  been  made  using  a  grid  of  113  x  30  points.  Such  a  grid  requires 
approximately  4  secs  per  time  step  with  a  potential  of  0.4  secs  per  time  step 
if  code  vector izat ion  were  performed.  The  question  of  the  run  time  per  case 
then  rests  on  the  number  of  time  steps  to  convergence. 

Although  the  cascade  Navier-Stokes  analysis  can  be  used  both  for 
time-dependent  and  steady-state  flow  situations,  and  although  good  results 
for  time-dependent  flows  have  been  obtained  with  an  airfoil  version  of  this 
code  (Ref.  28),  the  focus  of  the  present  invest igat ion  is  flows  whose  steady 
solution  is  sought.  In  such  a  case,  it  is  not  necessary  to  accurately  follow 
the  transient  motion  and  indeed  it  may  be  advantageous  not  to  follow  the 
transient  motion  accurately,  if  this  accelerates  convergence  to  steady 
state.  The  present  approach  is  based  upon  this  concept  and  utilizes  the 
matrix  conditioning  technique  of  Refs.  29  and  30  to  accelerate  convergence  to 
a  steady  state. 

Using  the  techniques  described  in  Refs.  29  and  30  has  allowed  cascade 
calculations  to  converge  very  rapidly  even  for  low  Mach  number  subsonic  flows 
which  in  general  can  be  difficult  to  converge.  A  typical  convergence  history 
is  presented  in  Fig.  7  which  presents  residual  versus  time  step  number  for 
the  subsonic  cascade  calculation  corresponding  to  the  experiment  of  Hobbs, 
et  al  (Ref.  14).  As  can  be  seen,  the  maximum  residual  in  the  flow  field 
drops  by  three  and  one-half  orders  of  magnitude  in  seventy  time  steps.  The 
reason  for  the  sudden  jump  at  time  sixty  will  be  discussed  shortly.  The 
maximum  residual  is  defined  hy  the  maximum  imbalance  of  any  equation  at  any 
point  when  the  time-derivative  terms  are  omitted.  At  seventy  time  steps  the 
solution  has  essentially  stopped  changing  and  as  will  be  shown,  the  surface 
pressure  calculated  is  in  very  good  agreement  with  that  measured.  At  4  secs 
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CRAY  CPU  time  per  time  step  the  present  calculation  can  be  considered 
converged  in  Less  than  300  CPU  system  secs. 

The  reason  for  the  sudden  jump  in  residual  at  time  step  60  is  due  to  the 
method  in  wh ich  the  calculation  was  run.  The  calculation  was  initiated  with 
artificial  dissipation  parameters,  dx  and  dy,  equal  to  0.5  (see  the 

*  section  on  artificial  dissipation).  At  time  step  sixty  this  was  redued  to 
0.05  except  in  the  immediate  vicinity  of  the  leading  edge  x/c  <  .02  where  it 

*  was  kept  at  0.5  due  to  marginal  resolution  in  this  region. 

A  comparison  between  the  calculated  and  measured  surface  pressure  for 
this  case  is  shown  in  Fig.  8.  As  can  be  seen,  the  calculated  and  measured 
values  are  in  excellent  agreement.  Calculated  velocity  and  pressure 
coefficient  contours  are  shown  in  Figs.  9-11,  and  a  velocity  vector  plot  is 
shown  in  Fig .  12 . 

The  second  Navier-Stokes  subtask  considered  under  the  present  effort 
concerns  code  input.  Under  the  present  effort  many  of  the  items  previously 
required  to  run  a  case  have  been  automat ical ly  set  thus  allowing  a  new  user 
to  run  the  code.  With  the  current  cascade  code  once  the  grid  has  been 
generated  only  the  following  items  are  required  to  run  the  code: 

1  -  Restart  flag 

2  -  Number  of  grid  points  in  each  direction 

3  -  Upstream  flow  angle,  upstream  total  pressure, 

downstream  static  pressure 

4  -  Reference  velocity,  temperature  density  and  viscosity 

to  set  reference  Reynolds  and  Mach  number 

5  -  Artificial  viscosity 

6  -  Number  of  time  steps  requested 
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ESTIMATES  OF  TECHNICAL  FEASIBILITY 


Under  the  present  effort,  an  existing  coordinate  generator /Navier-Stokes 
solver  has  been  automated  and  streamlined  so  as  to  allow  a  new  user  to  obtain 
accurate  cascade  solutions  in  an  efficient  manner.  The  coordinate  generator 
requires  relatively  little  input  other  than  blade  shape  and  spacing  input  and 
has  been  used  to  generate  two  very  different  cascade  coordinate  systems. 
Further  effort  of  this  type  for  the  *0*  grid  generator  would  lead  to  a 
robust,  easily  used  cascade  coordinate  generation  code  for  both  'O'  and  'C1 
grid  types . 

The  second  general  item  investigated  focused  upon  the  Navier-Stokes 
code.  The  work  done  under  the  present  effort  has  allowed  steady  cascade  flow 
fields  to  be  generated  in  less  than  300  secs  of  CRAY  CPU  time  without  resort¬ 
ing  to  code  vect or izat ion .  At  commercial  computer  rates  this  represents  a 
very  practical  cost.  Although  code  vector i zat ion  would  be  a  time-consuming 
procedure,  it  is  straight-forward  and  it  is  estimated  that  a  vectorized  code 
could  obtain  solutions  in  30  secs  of  CRAY  CPU  time.  This  vectorized  run  time 
would  reduce  costs  to  the  point  where  multiple  runs  on  a  daily  basis  could  be 
considered  if  so  desired  thus  allowing  the  Navier-Stokes  code  to  be  a  prac¬ 
tical  design  tool.  Extension  of  the  current  rapid  run  time  version  to 
include  an  energy  equation  would  be  a  straight-forward  task.  In  this  regard 
it  should  be  re-emphasized  that  the  current  code  does  include  an  energy 
equation  (although  not  optimized)  and  has  been  successfuly  used  in  heat 
transfer  calculations  (Ref.  16). 

The  present  effort  has  concentrated  upon  steady  flow  fields  at  very  low 
subsonic  through  transonic  Mach  numbers,  although  the  code  solves  the 
time-dependent  equations  and  could  be  used  for  time-dependent  problems. 
Results  obtained  for  flow  about  an  isolated  airfoil  oscillating  through  a 
dynamic  stall  loop  have  been  very  encouraging  in  verifying  the  time-accuracy 
of  the  procedure  for  free  stream  Mach  numbers  at  or  above  0.30  (Ref.  28).  It 
should  be  recalled  that  for  steady  flows  converged  solutions  can  be  obtained 
easily  for  much  lower  free  stream  Mach  numbers,  however,  the  technique  used 
is  not  time-accurate.  Current  efforts  at  SRA  are  aimed  at  modifying  the 
algorithm  to  allow  time-accurate  calculations  for  low  free  stream  Mach 
numbers . 
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Since  many  time-dependent  problems  require  small  physical  steps,  many 
time  steps  may  be  required  to  model  a  time-dependent  process.  In  these 
cases,  a  fast  vectorized  code  would  obviously  be  a  major  advantage.  Once 
such  a  code  was  available,  it  would  be  practical  to  apply  the  code  to 
problems  such  as  unsteady  flow  in  cascades  due  to  blade  wake  passing  and 
flows  with  back  pressure  or  inlet  perturbations  on  a  routine  basis.  The  code 
could  also  be  used  in  a  full  nonlinear  approach  to  the  rotating  stall 
problem. 


•t 
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Fig.  1  -  Constructive  coordinate  system. 


Fig.  2  -  Constructive  coordinate  system  -  outer  loop  parameterization. 
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Fig*  3  -  Constructive  coordinate  system  -  leading  edge  outer  loop  parameterization. 
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Fig.  4  -  Constructive  coordinate  system  -  secondary  loops. 


31 


Fig.  5  -  Coordinate  system  for  Hobbs  cascade. 
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Fig.  6  -  Coordinate  system  for  Turner  cascade. 
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Fig.  7  -  Iteration  history,  Hobbs  subsonic  cascade. 
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O,  □,  A  DATA  (REF.  14) 
-  CALCULATION 


Fig.  8  -  Calculated  and  measured  pressure  distribution, 
subsonic  compressor  cascade. 
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Fig.  9  -  U-velocity  contours  for  Hobbs  cascade. 
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Fig.  10  -  W-velocity  contours  for  Hobbs  cascade. 
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Fig.  11  -  Pressure  coefficient  contours  for  Hobbs  cascade. 
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Fig.  12  -  Velocity  vector  plot  for  Hobbs  cascade. 


APPENDIX  A  -  Splil  LBI  Algorithm 


Li  neariz.it  ion  and  Time  Di  f  f  orenc  i  ng 

The  system  of  governing  equations  can  be  written  as  a  single  grid  point 
in  the  following  form: 


3  H(<j>)/3t  =  D(<f>)  +  S(«f>)  (A-l) 

where  <f>  is  the  column-vector  of  dependent  variables,  H  and  S  are  column-vector 
algebraic  functions  of  ,  and  D  is  a  column  vector  whose  elements  are  the  spatial 
differential  operators  which  generate  all  spatial  derivatives  appearing  in  the 
governing  equation  associated  with  that  element. 

The  solution  procedure  is  based  on  the  following  two-level  implicit  time- 
difference  approximations  of  (A-l) : 

(Hn+1  -  Hn)/At  =  B(Dn+1  +  Sn+1)  +  (1-B)  (Dn  +  s")  (A-2) 

,  r  t  ,,n+l  ,  ,1/Jn+1N  ,  A  n+1  n  0 

where,  for  example,  H  denotes  H(<f>  )  and  At  =  t  -t  .  The  parameter  p 

(0.5  -  B  -  1)  permits  a  variable  time-centering  of  the  scheme,  with  a  truncation 
2 

error  of  order  [At  ,  (B  -  1/2)  At]. 

A  local  time  linearization  (Taylor  expansion  about  (f>n)  of  requisite  formal 
accuracy  is  introduced,  and  this  serves  to  define  a  linear  differential 
operator  L  such  that 

Dn+1  =  Dn  +  Ln  (4>n+1  -  <}>n)  +  0  (At2)  (A-3) 


Similarly , 


Hn+1  =  Hn  +  (3H/3<+On  Un+1  -  4>n)  +  0  (At2) 

Sn+1  =  s"  +  (3S/34>)n  (<J.n+1  -  <}")  +  0  (At2) 

Equations  (A-3, 4, 5,)  are  inserted  into  Eq .  (A-2)  to  obtain  the 
which  is  linear  in 


(A- 4) 

(A- 5) 

following  system 


(A  -  BAt  Ln)  (4>n+1  -  4>n)  =  At  (Dn  +  Sn)  (A-6) 

and  which  is  termed  a  linearized  block  implicit  (LBI)  scheme.  Here,  A  denotes 
a  square  matrix  defined  by 


A  =  (3H/3<Ji)n  -  BAt  (3S/3cf>)n 


(A- 7) 


Special  Trent mont  of  Diffusive  Terms 

The  time  differencing  of  diffusive  term?;  is  modified  to  nccomodnt  e 

cross-derivative  terms  and  nlso  turbulent  viscosity  and  artificial  dissipation 

coefficients  (i.e.,  p  ,  d.)  which  depend  on  the  solution  variables. 

e  1 

These  diffusive  coefficients  are  evaluated  explicitly  at  t  during  each 
time  step.  Nota t ionally ,  this  is  equivalent  to  neglecting  terms  proportional  to 
Dp ^ / S 4>  or  9d^/9<})  in  L™,  which  are  formally  present  in  the  Taylor  expansion  (A-3) , 
but  retaining  all  terms  proportional  to  or  d^.  in  both  L™  and  Dn. 

In  addition,  the  viscous  terms  in  the  present  formulation  include  a 
number  of  spatial  cross-derivative  terms  which  are  evaluated  explicitly  at  tn. 

To  preserve  notational  simplicity,  it  is  understood  that  all  cross-derivative 
terms  appearing  in  Ln  are  neglected  but  are  retained  in  Dn. 


Consistent  Splitting  of  the  LBI  Scheme 

To  obtain  an  efficient  algorithm,  the  linearized  system  (A-6)  is  split 
using  ADI  techniques.  To  obtain  the  split  scheme,  the  multidimensional 
operator  L  is  rewritten  as  the  sum  of  three  '’one-dimensional"  sub-operators 
(i  =  1,  2,  3)  each  of  which  contains  all  terms  having  derivatives  with 
respect  to  the  i-th  spatial  coordinate.  The  split  form  of  Eq .  (A-6)  can  be 
derived  either  by  following  the  procedure  described  by  Douglas  and  Gunn  in 
their  generalization  and  unification  of  scalar  ADI  schemes,  or  using  approximate 
factorization.  In  either  case,  for  the  present  system  of  equations  the  split 
algorithm  is  given  by 


(A  - 

BAtlp 

k 

(<f>  - 

4*n) 

=  At 

(Dn  + 

sn) 

(A-8a) 

(A  - 

BAtL") 

** 

(♦  - 

<pn) 

=  A 

■k 

u  - 

*n) 

(A-8b) 

(A  - 

BAtlp 

un+1- 

■<pn) 

=  A 

** 

(4-  - 

*n) 

(A-8c) 

k  k  k 

where  <p  and  <j>  are  consistent  intermediate  solutions  [18,20].  If  spatial 
derivatives  appearing  in  L_^  and  D  are  replaced  by  three-point  difference 
formulas,  then  each  step  in  Eqs.  (A-8a,  b,  c)  can  be  solved  by  a  block- 
tridiagonal  elimination. 
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